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PRECEDING PAGE BLANK NOT FILMED. 

FOREWORD 

A workshop meeting on numerical analysis was held at the National Aeronautics and Space 
Administration Headquarters, Washington, D.C., on January 29 and 30, 1968. Sixty persons, about 
half those invited as actively interested in numerical analysis applications to aerospace problems, regis- 
tered attendance during the 4 sessions of the meeting, and of these, 34 contributed papers of which 
abstracts are printed here. 

Nearly all of the contributed papers represent research done either at NASA centers, or elsewhere 
on contract or grant from NASA. Therefore, the spectrum of the collected abstracts shows many mathe- 
matical approaches to a wide range of aerospace problems, with perhaps the only essential characteristic 
being analysis directed toward ultimate numerical evaluation. 

Of the many, besides the authors, who contributed to the success of this meeting, C. V. L. Smith 
(Atomic Energy Commission, Germantown, Md.), J. W. Siry (Goddard Space Flight Center, NASA), and 
L. C. Young (University of Wisconsin), who acted as chairmen of the last three sessions, deserve special 
thanks. Finally, it should be noted that in planning and other managerial details, the active collaboration 
of R. N. Seitz and others at Marshall Space Flight Center was practically indispensable. . 

Raymond H. Wilson, Jr. 

Chief, Applied Mathematics Program, NASA 



NUMERICAL INTEGRATION OF LUNAR MOTION 
EMPLOYING A CONSISTENT SET OF CONSTANTS 


Charles J. Devine, Jr. 

Jet Propulsion Laboratory 

A numerical integration and least-squares fit to the current JPL lunar ephemeris LE 4 (contained 
in JPL development ephemeris DE 19) over a 2-year period from April 25, 1966, to April 26, 1968, is 
described. The method of integration used and the differential equations solved, along with numerical 
methods employed to reduce computational error, are discussed. The fitted lunar ephemeris, when 
compared with Lunar Orbiter ranging data, reduced the residuals by several hundred meters. 


REFERENCES 

Devine, C. J.: JPL Development Ephemeris No. 19. JPL TR 32-1 181, Nov. 15,1967. 

Mulholland, J. D.; Devine, C. J.; Gary, C. N.; and Sjogren, W. L.: Gravitational Inconsistency in the 
Lunar Theory. Part I— Numerical Determination. Part 1 1 — Confirmationvia Radio Tracking. 

To appear in Science. 



USE OF THE CONVOLUTION INTEGRAL 
FOR SIMULATION OF CONTINUOUS DYNAMICS 


H. H. Trauboth 

Marshall Space Flight Center, NASA 

A numerical method using the discrete convolution for the computation of all time responses of 
any complex multiloop control system is described in principle. The control systems are represented 
by block diagrams which are formed by the linear and nonlinear blocks well known in control engineer- 
ing. The properties of these blocks are given by simple relations between input and output signals in 
the time domain at discrete and equal intervals. Any number of input signals may be fed into the sys- 
tem at any point of the system. The input signals and characteristics of the blocks may be arbitrary 
time functions given in analytical or numerical form. 

For each linear block and each time interval the convolution must be calculated; this can be 
approximated by recursive formulas with a fixed number of only a few arithmetic operations. The 
results are then input to a system of algebraic equations which can be solved directly in the case of 
linear physical systems or iteratively in the general nonlinear case. 

The method has the advantage that the block diagram can be directly mapped into a matrix equa- 
tion where the position of the coefficients corresponds to the position of the blocks in the block dia- 
gram. There is no need for special treatment of algebraic loops and of sequencing. The mathematical 
method is especially suited for programing an electronic digital computer. 

The method proved very powerful and fast compared with conventional methods of solving ordi- 
nary differential equations in a number of applications. 

It is felt that the convolution method deserves more attention by numerical mathematicians in 
order to make optimum use of it and to explore its theoretical limitations. 


% 

REFERENCE 

Trauboth, H. H.: Digital Simulation of General Control Systems. Simulation, June 1967. 
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THE CALCULATION OF ATMOSPHERIC PROPERTIES 


Wm A. Mersman 
Ames Research Center, NASA 

Current methods of computing the properties of the Earth's atmosphere are quite efficient'for 
altitudes less than 90 kilometers. For higher altitudes the methods are intolerably slow, primarily 
because of the difficulty encountered in computing the pressure. This requires the integration of a 
rational fraction of high degree, based on the exact formula for geopotential altitude in terms of 
geometric altitude. 

A simple, linear rational fraction for geopotential altitude as a function of geometric altitude has 
been obtained by combining analytical and numerical techniques. 'Use of this approximate formula 
provides an algorithm for calculating the pressure with a speed that is essentially independent of alti- 
tude. Thus, above 90 kilometers, the new method is about 12 times as fast as current methods, with 
no loss in accuracy. 
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THE SYSTEM OF SINGULAR INTEGRAL AND ORDINARY DIFFERENTIAL EQUATIONS 
FOR THE COMBINED RADIATION CONDUCTION HEAT TRANSFER 
OF A FINNED SPACECRAFT CYLINDER 


Edward M. Keberle 
Jet Propulsion Laboratory 

The computation of the equilibrium temperature distribution function of a finned cylinder in the 
radiation field of the Sun takes into account the following energy transfers at each infinitesimal surface 
element: 

(1) Absorption and reflection of incident radiation from the Sun or other parts of the finned 
cylinder 

(2) Emittance of radiation from the surface according to the relationship W= eaT 4 

(3) Conduction of heat inside the surface to neighboring surface elements. 

These requirements lead to the followingsystem of equationsin terms of the unknown temperature 
T(x) as well as of the unknown radiosities B(x) resp. C(x) of the outer resp. inner parts of the cylinder: 

d!l£x) +^l ilM + a t 4 (x)+B(x)+ a C ( x )+ a = o 
dx 2 x dx 2 4 d 

B(x) + |3 1 T 4 (x) + /3 2 f ll+l2+l * K uter (x,y)B(y)dy+0 3 = 0 

C(x) + y t T 4 (x) + y 2 f 1 2 K. nner (x,y)C(y) dy + y 3 = 0 

* 7 / I 

l 


The mathematical complexities consist of - 

M 

(1) A combined system of nonlinear ordinary differential and nonlinear integral equations 

(2) Singular kernel 

(3) Nonlinear boundary conditions 

(4) Discontinuous solution B(x), C(x) 

Such a problem cannot be solved either analytically or by the construction ofparticular subroutines. 
Relevant is the applied mathematical and numerical task of adjusting common machine subroutines to 
the full realistic problem, which is filled with traps and exotic requirements. 

A variant of the method of successive approximations was used to construct B n (x),C n (x) by 
means of the integral equations while keeping T n _ x (x) fixed and then solving the boundary value prob- 
lem of the differential equation to obtain T n (x), while B n (x),C n (x) remain unchanged. 

Convergence was improved and accuracy of 0.1 percent obtained by internal repetition of the 
differential equation routine in comparison with the integral equation routine. 

REFERENCES 

Bunton, W. R.; and Keberle, E. M.: Space Programs Summary, 37-42, vol. IV. Supporting Research 
and Advanced Development, Jet Propulsion Laboratory. 

Plamondon, J. A.; and Horton, T. E.: Intern. J. Heat Mass Transfer, vol. 10, pp. 665-679, 1967. 
Sparrow, E. M.; and Cess, R. D.: Radiation Heat Transfer, 1966. 
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THE AUTOMATIC INTEGRATION OF 
STIFF ORDINARY DIFFERENTIAE EQUATIONS* 


C. W. Gear 

University of Illinois, Department of Computer Science 

Stiff differential equations arise in physical systems in which two or more time constants are greatly 
different. They can arise also when a single time constant is very small in relation to a forcing function. 
An example of the latter is the equation y'= X(y- F(t)) + F'(t), which has the solution y = ce^ t +F(t). If 
\«0 and F(t) is a smooth function, the quantity ce^ is negligible after a short time period. There- 
after, it would be desirable to use step sizes h that control the error due to F(t). Unfortunately, the 
stability of the numerical methods in current use requires that | hX | be bounded (typically it must be 
less than 1). We now have multistep methods that are stiffly stable in the sense that they are stable for 
all Re(hX)<D <0 for some D (ref. 1). These methods have allowed the step size to be increased by 
several orders of magnitude over the conventional methods such as Adams. 

A class of multistep methods are used for the integration (refs. 2 and 3). These references show 
estimates of the various derivatives available. The difference between the predictor and the corrector 
is also used to estimate higher derivatives. These estimates make it possible to change the step size and 
the order so that a minimization of the amount of work per step might be attempted. An additional 
advantage of automatic order selection is that the initial order can be set to one. Because it is then a 
one-step method, no special starting procedures are required. Comparisons of these methods with 
Adams’ methods have been made for nonstiff, moderately stiff, and very stiff equations (ref. 4). The 
numerical results indicate that Adams’method is better than the stiff method by a factor of less than 
Jl for nonstiff equations, but very much inferior for stiff equations. The automatic order methods are 
superior to the fixed order methods for the equations considered. 

It is concluded that stiff methods should probably be used if the question of stiffness is open. If 
fl the user can be certain that his equations are not stiff, an Adams-type method is better. One unfortu- 
nate feature of the method is that the partial derivatives dyjdy^ are needed approximately. To save 
the user the work of preparing subroutines for both y'(y) and 3y'(y)/ 3y, a package is being constructed 
which will accept the differential equations in symbolic form and compile 360 codes for the functions 
y ‘(y) and their partial derivatives. It is proposed to include the Adams’ and the stiff methods in the 
final package so that the choice of method, order, and step can be made. 

REFERENCES 

1. Numerical Integration of Stiff Ordinary Differential Equations. Rept. no. 221. Department of Com- 

puter Science, Univ. of Illinois, Urbana, 111. 

2. Gear: Report no. ANL-7 126, Argonne National Laboratory. 

3. The Numerical Integration of Ordinary Differential Equations. Math. Comp., vol. 21, no. 98, pp. 146- 

156. 

4. Automatic Numerical Integration of StiffOrdinaryDifferentialEquations. Submittedto IFIPS, 1968. 


*Work performed under grant ABC (1 1-1)- 1469 from the U.S. Atomic Energy Commission. 
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SERIES REPRESENTATIONS FOR GENERALIZED TEMPERATURES 


Deborah Tepper Haimo 
Southern Illinois University 

The generalized heat equation is given by 


A x u(x,t) = ^ u(x,t) (1) 

where A x f(x) = f”(x) and v is a fixed positive number. A generalized temperature is defined 

as a C 2 solution of equation (1). The fundamental solution of equation (1) is the function 
G(x;t) = G(x,0;t) where 


G(x,y;t) 


er( x2+ y 2 )/ 4t ^ 

(2t) u+1 / 2 2t 


with J2(z) = 2 u-1/2 r(i;+l/2)z 1 ^ 2_lJ I u _ 1/2 (z), I a (z) being the ordinary Bessel function of imaginary 
argument of order a. The solution of the Cauchy problem of equation (1) with u(x,0) = x 2n is the 
polynomial 


p n / x - t > = £ 2 

k=0 


2k A r(u+ 1/2 + n) 

lk/r(u + l/2 + n-k) 


x 2n-2k t k 


The Appell transform of P njU (x,t) is the function 


V 


W ni „(x,t) = G(x;t)P n> „^-A 

A function u(x,t) belongs to the class H* for a<t <b if and only if it is a generalized temperature 
satisfyingthe condition that for all t, t', a<t'<t<b, 


J s*QO _l 

! G(x,y;t-t)u(y,t') d/x(y), d/r(y) = [2 u ~ 1/2 r(u + l/2)] y 2u dy, 

0 

the integral converging absolutely for all x. Note that the polynomial P n y(x,t) belongs to H* for 
-co<t<co ; and its Appell transform W n y (x,t)eH* for t >0. 

Theorem: A necessary and sufficient condition that 


OO 


u(x.t) =2 a n P n ) i>( x ’ t ) 

n=0 
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for | t|<a is that u(x,t) eH* for | t|<a. The coefficients a n are given by either 


_ u^ 2n )(0,0) 
a n * (2n)! 


or 


a„ - 


n - (fan) I, W n ,„(*.t>«(x,-t)d»(x), 0<t<o. 


Theorem: A necessary and sufficient condition that 


OO 


u(x,t) = £ b n W n 

n=0 


for t >o >0 is that there exist a function v(x,t) eH* for i 1 1 < 1/a such that u(x,t) will be its Appell 
transform. Here b n = v^ 2n ^(0,0)/(2n)! 

The theory for the representation of generalized temperature functions in series of polynomial 
solutions and of their Appell transforms is analogous to that for the expansions of analytic functions 
in series of monomials p n (x) = x n and of inverse monomials w n = l/x n+1 , respectively. 


REFERENCES 

Haimo, D. T.: Expansions in Terms of Generalized Heat Polynomials and of Their Appell Transforms. 

J. MathMech., vol. 15, no. 5, May 1966, pp. 735-758. 

Haimo, D. T.; and Cholewinski, F. M.: Classical Analysis and the Generalized Heat Equation. To 
appear in SIAM Review. 



APPLICATION OF TWO TECHNIQUES IN THERMAL RADIATION PROBLEMS 


P. A. Newman and R. W. Barnwell 
Langley Research Center, NASA 

The two techniques which are applied to thermal radiation problems are Lax’s finite difference 
method and Case’s singular normal mode expansion method. Lax’s method is a numerical technique 
for solving nonlinear hyperbolic equations, whereas Case’s method is basically a rigorous mathematical 
technique for solving the linear transport equation. Neither of these techniques is new, and they are 
not applied to the same type of radiation problem. However, applications involving thermal radiation 
and emission are recent and are currently being investigated. 

Specifically, Lax’s finite difference technique is used to solve the problem of nonadiabatic flow of 
a thermally radiating gas past a sphere. The solution of the steady-state flow problem is sought as the 
asymptotic solution of the time-dependent flow problem; this method has been used in the past for 
solving steady-state flow problems in this way. However, when the gas is allowed to absorb and emit 
radiation, integral terms appear in the energy equation. It is shown that these can be treated in a 
straightforward manner, since the complete solution is known on time planes previous to the one on 
which the solution is being sought. Therefore, integrals across the flow field can be evaluated. 

The second technique was developed by Case to solve one-speed, one-dimensional neutron trans- 
port problems. Basically, he showed that a separation of variables in the homogeneous linear transport 
equation could be performed and that the resulting eigensolutions, although singular, were complete 
and orthogonal. Thus, any other solution can be expanded in terms of them, the expansion coefficients 
being determined by the source and boundary conditions of the particular problem. Application of 
Case’s technique in radiative transfer is illustrated by solving the problem of a finite slab of gray gas 
between infinite parallel plates. The use of this technique allows one to obtain accurate solutions with 
very little numerical analysis. 
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BAIRSTOWS METHOD REVISITED 


R. W. Hamming 

Bell Telephone Laboratories, Inc. 

The basic purpose of this report is to point out the desirability of designing library routines to 
meet the use that is going to be made of the answers, rather than merely meeting the obvious mathe- 
matical criteria. Bairstow’s method for finding zeros of polynomials is used as an example. 

Among other points made is that a routine which reports close distinct zeros will usually cause 
severe roundoff problems much later in the computation, whereas if they are reported in a multiple 
zero then the roundoff will not occur. 

The paper discusses methods of overcoming some of the well-known faults of Bairstow’s method 
and identification of multiple zeros. Simplicity in understanding what is going on is also emphasized. 
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ANALYTIC APPROXIMATIONS TO THE SOLUTIONS 
OF VARIATIONAL EQUATIONS 


Samuel Pines 

Analytical Mechanics Associates, Inc. 

The solutions of the differential variational equations of a dynamical system are often required 
in problems involving gradient search methods, two-point boundary value problems, optimal control 
theory, differential correction schemes, and so forth. Usually this is accomplished by numerical inte- 
gration at the cost of an increase in complexity and computing time. Since, for many search problems, 
the exact value of the partial derivative is not specifically required, a method is presented for obtaining 
an approximation to the solution of the variational problem by analytically differentiating a known 
closed-form solution of an approximating dynamical system. 

Specific examples are described, such as the closed-form solution of the two-body problem for 
perturbation problems in orbital mechanics and the known solution of the motion of a vehicle in a 
uniform gravitational field under the action of an inertially fixed thrust of constant magnitude for use 
in error propagation of vehicles under thrust. 
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QUASILINEAR FORNULATION OF NONLINEAR TRANSFORMATIONS 


Henry Wolf 

Analytical Mechanics Associates, Inc. 

In many computations, a set of parameters has to be changed. Generally, two types of transfor- 
mation formulas are used— one for infinitesimal and one for finite changes. The infinitesimal trans- 
formation is usually linearized and behaves properly when the changes tend to zero, but behavior 
becomes improper when the changes become finite. The finite transformations, on the other hand, do 
not properly tend to zero when the transformation is small. The boundary between the two regimes 
where each transformation is valid is usually ill defined. Further, transformations often depend on 
several parameters, some infinitesimal and some finite. In this case, neither formula is completely 
satisfactory. 

In the paper a method is presented to obtain a “quasilinear” form, which is valid for finite changes, 
but contains the nonlinearities as a coefficient of the transformation parameters, and consequently 
behaves properly also for the infinitesimal case. 

Applications of this method are presented. They include quadratic Newton-Raphson iteration, 
the planetary and Encke perturbations in equations of motion, vector rotation formulas, and differen- 
tial correction schemes. 
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IS THERE A UNIFIED SOLUTION FOR A CUBIC EQUATION? 

Mary H. Payne 

Electronics Research Center, NASA 

A method of quickly and accurately solving the cubic equation ax 3 +bx 2 +cx+ d = 0 is desired for 
various combinations of “largeness” and “smallness” of the coefficients. If, for example, a is not too 
small, the standard Cardan method is acceptable. If a is small but d is not, one can divide by x 3 and 
solve for 1/x. If 3ax 2 + 2bx + c is not too small, Newton’s method can be used. The question posed is, 
How does one formulate a sequence of tests on the coefficients together with methods of solution for 
various outcomes of these tests so as to cover all cases? Suggestions made by Gear (University of Illi- 
nois) and Hamming (Bell Telephone Laboratories) are currently under study. 
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THE NUMERICAL SOLUTION OF SINGULAR INTEGRAL EQUATIONS 

OF POTENTIAL THEORY 

M. S. Lynn and W. P. Timlake 
IBM Houston Scientific Center 

The authors discuss integral equations that arise out of three-dimensional potential theory; i.e., 
where the kernel is* 

Equations of this type are singular both in the sense of the kernel being singular [a( 1/r) ] as well 
as the operator being singular. Thelatter singularity is deflated out by seeking that particular solution 
<p for which 

J pdA = 0. 

The surface S is “triangulated” into n areal elements of maximum diameter p . A numerical scheme 
derives a system of approximating linear algebraic equations of order n in ^ . 

Based on the assumption, inter alia, that the surface S is regular in the sense of Kellogg and 
Lyapunov, it is shown that \p differs from the average value of <p in the norm by at most Kp^ +1 ^ 3 
(K is a constant). The algebraic equations are solved by Jacobi iteration, in which the spectral radius is 
bounded (as n” 00 ) strictly away from 1 ; for example, for the model problem, a sphere, the spectral 
radius is of order 1/2. Extensions to multiconnect regions were also considered. 


*For notation, precise statements, and a proper bibliography, see a paper with the same title which is to appear 
shortly inNumerische Mathematik. 
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SOLUTION OF SPARSE LINEAR SYSTEM OF EQUATIONS 

R. P. Tewarson 
State University of New York 

The problem of transforming a given sparse matrix A into a block diagonal form (bdf) and the 
subsequent transformation of each of the block diagonal matrices into as nearly upper triangular form 
(utf) as possible, by using only row and column permutations, is discussed. It is shown how some of 
the results from graph theory can be used to transform A to the bdf. To transform the block diagonal 
matrices into the utf’s, two methods are described, one of which uses linear programing while the 
other uses approximate probabilistic arguments. The latter method, in relation to the computational 
effort, yields significant results in practice. 
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AUTOMATIC ERROR BOUNDS FOR REAU ROOTS OF POLYNOMIALS 
HAVING INTERVAL COEFFICIENTS 


Richard J. Hanson 
Jet Propulsion Laboratory 

The notion of interval arithmetic (ref. 1 ) is briefly described, and a Newton-Raphson algorithm 
using interval arithmetic and range extensions of functions is discussed (ref. 2). This algorithm is 
particularly useful for fiiding real roots of polynomials with interval coefficients, and several examples 
are presented. 

The discussion is concluded with a presentation, and example, of an algorithm for solving a square 
linear system of algebraic equations of full rank (ref. 3). 


REFERENCES 

1. Moore, R. E.: Interval Analysis. Prentice-Hall, 1966. 

2. Hanson, R. J.: Automatic Error Bounds for Real Roots of Polynomials Having Interval Coefficients. 

To be published. 

3. Hansen, E.; and Smith, R.: Interval Arithmetic in Matrix Computations, part II, J. SIAM Num. 

Analysis, series B, vol. 4, no. 1, 1967, pp. 1-9. 
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NONLINEAR PROBLEMS AND 
MATHEMATICAL METHODS IN AEROELASTICITY 


E. H. Dowell 
Princeton University 

A discussion of methods for mathematically treating the dynamic interaction of a flowing fluid 
and flexible solid (aeroelasticity) is given. Primary emphasis is placed upon techniques for handling the 
nonlinear behavior of such physical systems. The choice of numerical procedures is based upon the 
great ability of high-speed digital computers to numerically integrate systems of ordinary integral- 
differential equations (initial value problem). Hence, the nonlinear partial differential field equations 
governing the fluid and solid motions are reduced to such systems by utilizing generalized Fourier or - 
eigenfunction expansions associated with the linear dynamics of the solid. Given these expansions, the 
fluid motion is determined from the fluid field equations, through the methods of the transform cal- 
culus, in terms of the (as yet unknown) time history of the solid motion. The fluid forces on the solid 
are then calculated, and the system of ordinary (nonlinear integral- differential) equations for the solid 
motion are subsequently determined from the solid field equations. Numerical integration of these 
ordinary integral-differential equations gives the time history of the solid motion. The flow field time 
history also may be determined if desired. Finally, it should be mentioned that the Fourier or eigen- 
function expansions for the linear solid which are needed may be determined analytically for simple 
solid geometries; however, for irregular geometries, the eigenfunctions must be calculated by yet an- 
other numerical procedure, e.g., the finite difference method. Less frequently, this also may be true of 
the flow field. It is emphasized that finite difference methods applied directly to the full nonlinear, 
fluid-solid interaction problem are impractical from the point of view of computer time and storage 
requirements. 
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REGULARIZATION AND STABILIZATION OF NUMERICAL SOLUTIONS 
OL DILLERENTIAL EQUATIONS 


Victor Szebehely 
Yale University 

The transformations which regularize the equations of motion of celestial mechanics involve the 
independent and the dependent variables. The purpose of these transformations is to eliminate the 
singularities occurring at collisions of two bodies; however, it is well known that even when collisions 
do not occur, the appearance of close approaches makes regularization mandatory from the point of 
view of numerical analysis. It was found that these regularizing transformations offer improved stability 
characteristics of the equations of celestial mechanics. 

Examples regarding the restricted problem of three bodies as applied to Earth-to-Moon trajectories 
show a significant time saving on electronic computers when regularizing variables are used. In com- 
parison studies, the methods of regularization show advantages over methods which do not eliminate 
the singularities in accuracy and in efficiency. Analytical approaches attempting to arrive at series 
solutions to Apollo-type trajectories also benefit from the process of regularization. 

Another example is the problem of integrating systems of ordinary differential equations of high 
order, say several hundred. Such n-body integrations show sensitivity and sudden increases in the 
error at any time when close approaches of the participating bodies occur. 

The final example demonstrating the importance of regularizing transformations is the solution of 
a problem of three bodies proposed by Burrau in 1913 and solved only recently. The repeated close 
approaches of the three participating bodies necessitate the use of regularizing variables originally pro- 
posed by Levi-Civita and applied to Burrau’ s problem by Szebehely and Peters. 


REFERENCE 
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Szebehely and Peters: Astron. J., vol. 72, Sept. 1967, pp. 876-883. 
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SUBARC ELIMINATION AND OPTIMAL ORBIT TRANSLER 


Thomas L. Vincent 
University of Arizona 

A procedure for handling optimization problems with corners resulting from the imposition of 
restrictions is presented. The requirements which must be met in order to apply the methods presented 
are that the dynamical equations of constraint contain only one control variable and that these equa- 
tions can be analytically integrated along any segment of the trajectory in which the control law is 
restricted. Under such circumstances it is shown that the restricted segments of the optimal trajectory 
may be effectively eliminated in the process of determining the changes in the Lagrange multipliers, 
state variables, and Hamiltonian function from one side of the restricted segment to the other. A set 
of conditions which may be used to determine the changes or “jumps” in the Lagrange multipliers and 
Hamiltonian function are obtained directly from a set of corner conditions developed for trajectories 
with discontinuities in the state variables. With the use of these conditions, there is no need to inte- 
grate the Euler equations along the restricted arc, as is usually done in the literature. 

As an example as to where these methods may be directly applied, the thrust-coast-thrust optimal 
orbit transfer problem is cited. It is shown that the optimal cutoff and relight points for the thrust 
vector could be obtained directly from the corner conditions, rather than by integration of the Euler 
equations along the restricted arc. 


REFERENCE 
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A SEGMENTAL EXTENSION OF SOLUTIONS TO OBTAIN A NUMERICALLY STABLE 

BOUNDARY VALUE PROBLEM SOLUTION 


J. P. Avery 

Arizona State University 

A technique is presented to circumvent the possible loss in accuracy associated with the super- 
position of initial value-problem solutions to obtain the solution to a specified boundary value problem. 
Should the dependent functions Y. of the initial value problems exhibit a rapid growth with increasing 
independent variable x, then serious cancellation of leading digits could result upon superposition to 
meet the remote edge boundary conditions. 

The suggested technique entails a repeated use of superposition of initial value problems at small 
segmental intervals in x in order to suppress the offending growth in the functions Specifically, at 
each segmental station x ; , initial value problem solutions are combined, first, to meet the desired re- 
mote edge boundary conditions artificially at %, and again to provide the required independent com- 
plementary solutions for continuing through the next segment. The resulting solutions then provide 
sets of initial conditions for extending the integration over the next segment to x^ , where the super- 
position procedures are repeated. The process is iterated until, at the remote edge x n , the boundary 
conditions are satisfied through the final segmental superposition of solutions. By keeping the segment 
size suitably small, the growth in magnitude of the dependent functions is accordingly controlled, and 
acceptable error is introduced at each superposition operation. 
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VARIABLE-MESH, MULTISTEP METHODS FOR 
ORDINARY DIFFERENTIAL EQUATIONS’’ 

R. Van Wyk 

North American Rockwell Corp., Rocketdyne Division 

A theory for variable-mesh, multistep, predictor-corrector solution of ordinary differential equa- 
tions is presented. The difference equations employed are generalizations, for the case of variable mesh 
spacing, of previous formulas requiring fixed step size. In addition to retaining the high local accuracy 
and convergence properties of the earlier methods, the variable mesh method is developed in a form 
conducive to the generation of effective criteria for the selection of subsequent step sizes in the step- 
by-step solution of differential equations. These criteria are based on considerations of truncation 
error, convergence of corrector iterations, and an extensive treatment of relative numerical stability. 
This treatment included tracking the roots of a cubic equation, whose coefficients are complex in the 
case of simultaneous differential equations, and determining stability as a function of one real and one 
complex parameter. The algorithm has been tested extensively and compared with other methods. 

The results of the comparison favor the new method. 

Current research on the extension of the variable mesh method to systems whose components 
decay at greatly different rates is also presented. 


: This work was supported by the National Aeronautics and Space Administration, Office of Advanced Research 
and Technology under Contract NAS7-470. Much of the material presented was taken from Rocketdyne Research 
Report No. 67-12, Sept. 1967. 
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STEP-SIZE SELECTORS FOR THE SHANKS FORMULAS OF RUNGE-KUTTA TYPE 


E. Baylis Shanks 
Vanderbilt University 

In 1962 and later, the author* developed the only Runge-Kutta-type formulas in existence of orders 
7, 8, 9, and 10. In addition, more efficient formulas were derived using 4, 5, 6, 7, and 10 evaluations. 
The 9th- and lOth-order formulas were developed in 1967. To indicate the complexity involved, it is 
noted that the lOth-order formula is derived by solving considerably more than 500 algebraic equations, 
more than half of which are of the 18th degree. 

A Computation Division report from Marshall Space Flight Center, NASA, states that Shanks’ 
formulas through the 8 th order require only about one-third as much machine computation in orbital 
and interplanetary probes as the previous Runge-Kutta-type formulas. 

Recently, the author has discovered that a second formula, less accurate but using no new evalua- 
tions, can be developed for each of the Shanks formulas. Let us indicate the approximations to y by 
Y x and Y 2 > the errors by e t and e 2 , and the evaluations of f by so that y 0 + h^c i f i =Y 1 = y+e 1 , 
where y 0 is the initial value of y and cq are constants of the formula. Subtracting Y 2 , we get 
h^( c i- dj) ft=Yi -Y 2 =ei-e 2 . The quantity S = | h^Oq-dj) fj | is the absolute value of a fairly good 
approximation to the error e x of the original formula. It is called the step-size selector and can be 
used to change the step size for the next step of the integration process (if necessary) so that 

icr n <s<io - n+4 

for an appropriate integer n. Other intervals for S may be used. Tests indicate that, in order to main- 
tain the same accuracy, formulas using the selectors will require only about one-third the computing 
time used for the original formula alone. 

These innovations are programed and in use at the Computation Division of Marshall Space Flight 
Center. 


REFERENCES 

Shanks, E. B.: Formulas for Obtaining Solutions of Differential Equations by Evaluations of Functions, 
Department of Mathematics, Vanderbilt University, Nashville, Tenn., 1963, pp. 1-12. 
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RATIONAL APPROXIMATION OF SPECIAL FUNCTIONS 


Henry C. Thacher,Jr. 

University of Notre Dame 

W. J. Cody 

Argonne National Laboratory 

Subroutines based on the best (minimum- maximum relative error) rational approximations are 
ordinarily the most efficient for computing mathematical functions. For continuous functions of a 
single real variable, several Fortran and Algol programs are available which can produce rational func- 
tions differing from the best approximation by no more than 0. 1 percent in less than 1 minute on even 
a moderate- speed computer. The major requirement for using these programs is a reference routine 
producing values of the desired function to three to four digits more than will be required for the 
approximation. Typically, the resulting approximations yield up to 8 significant digits for the ratio of 
2 cubics, and up to 20 for the ratio of 2 sextics. 

Aside from constructing and checking the reference routine, the principal difficulties are the 
choice of the particular auxiliary function to be approximated for various portions of the domain, and 
the determination of the intervals over which each form is to be used. 

J. F. Hart et al. show the best rational approximations ranging up to 25S for the elementary, and 
several higher functions, including the error function, the gamma function, the complete elliptic inte- 
grals, and the Bessel functions J Q (x) and JjCx). It also gives an extensive bibliography of other pub- 
lished approximations. 

Several other groups are also active in approximating additional functions. We would be interested 
in learning of other functions which are of sufficient importance in applications to justify their approx- 
imation. 


REFERENCE 
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EXTENSIONS AND APPLICATIONS OF THE HOUSEHOLDER ALGORITHM 
FOR SOLVING LINEAR LEAST-SQUARES PROBLEMS 


Charles L. Lawson 
Jet Propulsion Laboratory 

Let A denote a real matrix of rank r having m rows and n columns and let b denote a real 
vector having m elements. The problem-of determining an m-vector & which exactly or approxi- 
mately satisfies Ax =S occurs in many contexts. Wfe distinguish three cases: (l)m=n, (2)m>n, 
and (3) m<n, and further subdivide each case into two subcases: a, r=min(m,n) and b, r<min(m,n). 

In case ( 1 )a the matrix A is square and nonsingular. In the five other cases it is necessary to 
produce a square r x r matrix of rank r at some stage of the computation and solve a system of equa- 
tions involving that matrix. For example, the commonly used method of solving case (2) a in the sense 
of least squares involves the formation of the rxr matrix A T A. Similarly, a standard method for 
solving case (3)a uses the rxr matrix AA T . 

It is known that the spectral condition number of A T A (or of A A T ) is the square of the condition 
number of A. It is also known that there are methods using orthogonal transformations which permit 
the reduction of A to an r x r matrix whose condition number is the same as that of A. Because of 
the conservation of the condition number, these latter methods are, in general, capable of preserving 
more accuracy in the solution than the methods using A T A or A A T , assuming the same type (i.e., pre- 
cision) of arithmetic is used. More specifically, a computation which can be carried out satisfactorily 
with orthogonal methods using s-digit arithmetic may require 2s-digit arithmetic, if done using A A 
or AA t . 

In some, possibly most, cases of badly conditioned systems of linear equations, it is more meaning- 
ful to compute a solution of an appropriate neighboring system of lower rank that is reasonably well 
conditioned than to compute an accurate solution for the badly conditioned problem. Certain ortho- 
gonal methods adapt rather naturally to this approach. 

Computer users will be better served as subroutines using orthogonal methods become more 
readily available and these methods are more widely understood by people who specify and design 
computer programs. For example, reference 1 illustrates how a user's time can be absorbed in an 
essentially random search for a reliable least-squares program. Fortunately, one of the programs this 
user encountered happened to use orthogonal methods and, indeed, he did find that it was the only 
satisfactory program for his test cases. 

One algorithm for solving problems ( 1 )a and (2 )a by orthogonal methods, specifically by the use 
of Householder orthogonal transformations, is given in reference 2. In an effort to improve the flexi- 
bility and accessibility of orthogonal methods for general applications, we have extended the algorithm 
of reference 2 to provide a uniform approach to problems (l)a, (2) a, and (3 )a. With a somewhat arbi- 
trary criterion for rank estimation, this algorithm also encompasses problems (l)b, (2 )b, and (3 )b. 

This algorithm permits sequential processing of the rows of the augmented matrix [A :b] so that the 
storage required depends only on n and not on m. This algorithm has been implemented as a Fortran 
IV subroutine. 

A different use of orthogonal methods is the singular value decomposition algorithm of reference 
3. Wfe have converted the Algol code given in reference 3 to Fortran IV. This provides another 
numerically sound way to obtain the capabilities mentioned in the preceding paragraph. 
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It is expected that the ability to produce a numerically sound solution for Ax = ft regardless of 
the relationship between m, n, and r rail be valuable in many applications. We are investigating the 
use of these methods in the solution of nonlinear equations, nonlinear least squares, linear and non- 
linear Ip minimization (1 computation of independent eigenvectors associated with multiple 

eigenvalues, root locus computations, implicit function evaluation, minimization subject.to constraints, 
and solution of large systems of equations. 
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TRUNCATION ERRORS ON THE COMPUTATION OF ELLIPTIC INTEGRALS OF THE 
FIRST AND SECOND KINDS BY GAUSS AND LANDEN TRANSFORMATIONS 


Edward W. Ng 
Jet Propulsion Laboratory 

It is well known that algorithms based on quadratic transformations yield a high rate of conver- 
gence in the computation of elliptic integrals, which occur frequently in astronomical and astronautical 
problems, such as in artificial satellite theory. Four algorithms based on descending and ascending 
Gauss and Landen transformations have been investigated quite thoroughly, either in the trigonometric 
or algebraic versions, the latter commonly known as the AGM scale (ref. 1). 

In actual computer application, these algorithms present themselves in some form of an infinite 
iterative process, to be truncated after a finite number of steps. Some bounds and estimates obtained 
for errors involved in such truncation are reported. Such error analysis is to serve a twofold purpose: 
first, as a theoretical guide before computer programing, and second, to give actual numerical value of 
the errors during computer programing. A small tolerance parameter is chosen consistent with both 
the algebraic and trigonometric versions and applicable in all four algorithms. Error bounds in terms 
of this parameter are given for the computation of elliptic integrals of the first and second kinds. 
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APPLICATION OF INTEGRAL EQUATION METHOD TO 
NUMERICAL ANALYSIS OF BIHARMONIC PROBLEMS 


Kwan Rim 

University of Iowa, Department of Mechanics and Hydraulics 

This presentation is based on reference 1 and on subsequent work which is to be published shortly. 
The method of analysis is based on a general formulation of a biharmonic stress function in terms of 
two harmonic functions. The harmonic functions are each represented by single-layer potentials, each 
containing an unknown source density function. Then the boundary conditions yield two coupled 
integral equations closely related to those of the Fredholm type. Numerical solutions for the unknown 
source densities are obtained by assuming that the source densities are piecewise constant around the 
boundary of the domain, thereby replacing the integral equations by sets of simultaneous algebraic 
equations which can be effectively solved by using a high-speed digital computer. Once the source 
density functions are determined, stress in the elastic domain may be found by computing second 
derivatives of the stress function. 

Numerical results obtained by the author, as well as other investigators, indicate that the integral 
equation method is extremely accurate at a small distance away from the boundary, but that, in gen- 
eral, the accuracy of the numerical solution decreases as the point at which computations are performed 
approaches the boundary. Improvement of the accuracy of numerical results on and near the boundary 
is therefore identified as the central problem of the integral equation method. 

In the previous report (NASA CR-779), it was first shown that one of the three quantities to be 
computed at the boundary, namely, the total stress, can be evaluated directly from the approximate 
stress function. The research reported herein shows that greater accuracy in computed boundary stress 
components can be obtained from the numerical solution if the numerically determined, piecewise 
constant source density functions are replaced by functions which, with their first derivatives, are con- 
tinuous on the boundary. In contrast to other reported results, this modification permits one to com- 
pute unknown boundary stress components directly from the approximate stress function without 
recourse to such indirect techniques as extrapolation, numerical differentiation, and so forth. Hence, 
a meaningful comparison between the specified boundary stresses and those boundary stresses con- 
sistent with the approximate biharmonic function is now made possible. 

Numerous methods have been suggested for the solution of biharmonic boundary value problems, 
but few numerical results have appeared in the literature. Even where some results have been published, 
the method employed has usually been of limited application, especially dictated by the boundary con- 
ditions and geometry. Further improvement in the integral equation method, with its simplicity and 
generality, is likely to provide the most powerful means of attacking the biharmonic problems, thereby 
opening new avenues of advancement in such vital areas as elastoplasticity, viscoelasticity, second-order 
elasticity, and optimization. 
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MULTIDIMENSIONAL INTEGRALS 


R. Cranley 

Southwest Center far Advanced Studies 

A critical examination is made of a class of formulas which have been used for multidimensional 
integration. These formulas integrate all multivariate polynomials up to a fixed degree exactly, and 
the evaluation points are chosen in subspaces of low dimension in order to minimize their number. 

It is shown that when the evaluation points are chosen in this way, serious errors may be intro- 
duced when the formulas are applied to general integrands. 

The repeated Gauss formulas are interpreted as formulas of similar type. It is shown, using a multi- 
variate Taylor expansion, that the Gauss formulas handle many more terms of the expansion than a 
multivariate formula (using an equivalent number of points) and that these terms are generally of 
greater importance. 

It is concluded that repeated Gauss quadrature is still the best method for evaluating multi- 
dimensional integrals. 
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USES OF THE SUMMATION BY PARTS FORMULA IN COMPUTATION 


Tomlinson Fort 
Emory University 

The summation by parts formula, probably due to Abel, has been used throughout the years in 
various ways. It is shown in this report how it can be used in many places to abbreviate computation. 
Illustrations are given for the computation of 1 n 2, it, and e. Some of its strictly mathematical uses 
are pointed out; for example, a very simple proof of Euler’s theorem on the convergence of infinite 
series. This theorem is useful in computation. 



EVALUATION OF THE GEOPOTENTIAL FORCES FOR AN ORBITAL EPHEMERIS 

WITH GIVEN HIGH-ORDER HARMONICS 

F. J. Lerch 

Goddard Space Flight Center, NASA 

The number of terms in the geopotential forces are rapidly becoming quite large, and care must 
be taken in their evaluation on the high-speed computers for the numerical integration of an orbital 
ephemeris. The geopotential, normally given in spherical coordinates, will suffice for discussion since 
it contains the basic functions and form as the geopotential forces. This potential is described in terms 
of zonal, tesseral, and sectoral surface harmonic terms in which the harmonic coefficients are primarily 
derived from satellite data. Each harmonic term consists of separable factors in spherical coordinates 
of geocentric latitude, longitude, and distance, including an associated Legendre polynomial referred to 
by degree n and order m. The individual terms are double summed over n=0 to N and m=0 to n 
yielding (N+ 1) (N + 2)/2 total harmonic terms for the potential. It is estimated that for 1 -meter 
satellite accuracy, a complete set of coefficients of degree N = 18, yielding 1 90 terms, is required. Vari- 
ation of terrestrial gravity measurements averaged over 5° square sectors of the Earth tends to support 
the above 1 -meter figure. Interconrparisons of geodetic tracking systems, including the laser systems 
on the NASA GEOS-I satellite, are indicative of the need of improved gravity fields. An interesting 
figure, obtained by associating an individual coefficient with each 5° sector ( 100 000 sq. mi.) of the 
Earth, would require a complete potential of about N=40, or 800 harmonic terms, as there are two 
coefficients per term. The Smithsonian Astrophysical Observatory has recently derived a gravity field 
from combined satellite and terrestrial gravity measurements that has a complete set of harmonic co- 
efficients up through N = 15. 

A method of evaluating the gravity forces for computer computation is presented and employs 
recursive formulas in spherical coordinates for accumulating the (N+ l)(N+2)/2 harmonic terms. This 
is compared with a method using similar recursive formulas expressed in Cartesian coordinates, which 
is shown to be slightly less efficient in computer execution time. This is primarily because the potential 
harmonic terms are not expressiblein separable factors in the Cartesian coordinates. Another method 
is presented where the associated Legendre polynomial in each harmonic term of the above double sum 
is expanded in its individual terms in the Cartesian coordinates and combined with the longitudinal 
factor expanded in Cartesian coordinates; this is then represented by a double sum over these combined 
individual terms. This representation of the potential yields a fourfold sum with an extensive kernel. 

It is presented, since it has been developed for use, and it demonstrates that care must be taken in 
designing computational schemes because it requires exorbitantly large computer time on the IBM 360. 
With the use of the spherical recursive formulas in evaluating the gravity forces and with the use of an 
llth-order Cowell integration, it takes about 1 minute of UNIVAC 1108 computer time to generate a 
GEOST definitive orbital ephemeris over 1 day, and it is estimated that it will require less than this 
amount of time for a weekly orbital ephemeris on the IBM 360/95. The gravity field employed in 
these latter calculations and used in the GSFC orbital intercomparisons of the GEOS-I geodetic tracking 
systems is of degree N = 15. 
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CERTIFICATION OF LIBRARY SUBROUTINES 


L. Richard Turner 
Lewis Research Center, NASA 

One of the important questions of numerical analysis is the quality of the subroutines for the 
approximation of elementary functions and other service chores such as exponentiation provided in a 
higher level language package; for example, Fortran, provided by manufacturer XYZ or by your own 
shop. 

An analysis of the IBM System 360 library was made at the Argonne National Laboratory and 
published in their technical report, ANL7321, 

Using slightly modified techniques, much of this work has been repeated and also carried out for 
much of the IBM 7094 Fortran version 13. The Lewis Research Center program, written in IBM For- 
tran IV, produces a histogram of error distribution, identifies the arguments of extreme error in each 
test interval, and computes the root-mean-square value of the relative error. The selection of sub- 
intervals and generation of test data is reasonably flexible. A source of higher precision values for 
each argument is required. 


4 * 
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CLOSED-FORM SOLUTIONS FOR COAST-ARC TRANSITION MATRICES 
ASSOCIATED Will I OPTIMIZATION PROBLEMS 


Hugo L. Ingram 

Marshall Space Flight Center, NASA 

A Cartesian coordinate, closed-form solution of the two-body problem is available from the work 
of W. H. Goodyear. This solution makes- use of the regularizing variable ^ ,such that dt/d \p =R, and 
thus the singularities which are usually encountered at eccentricities of 0 and 1 are eliminated from 
the problem. Also, Goodyear develops a solution for the partial derivative matrix which indicates the 
effects on the coordinates of the position and velocity vector at a particular time with respect to changes 
in the coordinates of the initial position and velocity vector. This work has been extended to obtain a 
closed-form solution for the Lagrangian multipliers or adjoint variables across a coast arc and their 
associated transition matrices. This information can be used to construct a Newton's iteration to satisfy 
the boundary conditions for optimization problems and to construct perturbative numerical integration 
schemes for optimization problems. 



A VARIATIONAL BOUNDARY METHOD FOR STRESS ANALYSIS 


Z. M. Elias 

Massachusetts Institute of Technology 

A direct variational method for solving boundary value problems in stress analysis when a repre- 
sentation of the general solution of the field equations is known can be based on Reissner’s variational 
equation. The method involves numerical integration over the boundary of the domain. For a problem 
in classical three-dimensional elasticity, the functional to be made stationary is 

1 p? ) Ui<is -f£$ - u °h ° i,is 

For computing I, the Papkovich-Neuber representation of the displacement vector, u i? in terms of four 

harmonic functions <b. can be used. This is 
*1 

2Guj= (3-4u)0 i -x j 0 j i - 0 o i 


For the problem of plane elasticity, the functional is similar to I. For computation, the general solution 
for the displacement vector in terms of two analytic functions f(z) and g(z) can be used. This is 


Uj + iu 


2 



1 + U 

4 


zf'+ g 


In plate bending, the functional takes the form 



and can be expressed in terms of the deflection w whose general solution in terms of 2 analytic func- 
tions f(z) and h(z) has the form 

w= - ~ R e (zf - h) 

4D l-o 


This boundary method takes advantage of the availability of the general solution of the field equations 
and should require less labor than interior methods, such as finite difference and finite element methods. 
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SOME USES AND POSSIBLE USES OF FRECHET AND 
OTHER DERIVATIVES IN A BANACH SPACE 


H. Melvin Lieberstein 
Director of Research 
The Mathematics and Biology Corp. 

Wfe have tried to collate a number of ideas into a cogent program. Our first acquaintance with 
the use of Frechet derivatives in numerical analysis was in the G. Fichera method of finding mini- 
mum L 2 error bounds in approximating solutions of boundary value problems by the least-squares 
continuous (i.e., not finite difference) method (refs. 1 and 2). The next was the L. V. Kantorovich 
theorem on convergence of the Newton method in a Banach space. This theorem and a modification 
of it by M, Altman were used by a student, Alan Elcrat, to prove existence (Math. Biosciences, to 
appear) of regular solutions of a boundary value problem related to blood flow as I had formulated it. 
Specification of the initial velocity profile was avoided by replacing uniqueness with mean-square 
asymptotic uniqueness (giving uniqueness of a rapidly acquired steady state) with the hope that this 
would indeed make the proof of existence of regular solutions feasible. Dr. Elcrat finds a sequence of 
linear boundary value problems whose solutions converge to a solution of our nonlinear boundary value 
problem. By mapping from c 2 + a to C a (Holder continuous functions with a Leray norm), he forces 
certain differential operators to be continuous, and he is then able to utilize the results of Avner Fried- 
man's book on parabolic type to obtain existence for the linear problems. Wfe now plan to approximate 
solutions of the linear problems using least squares, accumulating error bounds of these and the relevant 
Newton method to obtain an overall bound for the error, but we intend to try for the fastest possible 
convergence properties before we start computing. Toward this end, we are examining asymptotic 
stability of the analogs 

dx/dt = - [P'(x)]' 1 P(x) or dx/dt = - [P , (x 0 )]' 1 P(x) 

of the Newton method and modified Newton method, where P: Ac X-»Y and X,Y are Banach spaces 
and P'(x) is the Frechet derivative of P. We may possibly integrate such equations numerically using 
Runge-Kutta in order to find an iterative method which is faster than the Newton method. Work in a 
forthcoming textbook (Harper & Row) on 

dx/dt = -(I+B _1 C)x+ B _1 b , 
an analog for the general linear iterative methods 

x k+ 1 = -B' 1 Cx k + B' 1 b , 


for linear systems 


(B+ C)x =Ax = b , 

which is gotten from it by a crude Euler first-order finite difference replacement, is encouraging in this 
direction. Here the analog is asymptotically stable if the iterative method converges, and its asymptotes 
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are solutions of the linear system. The converse is not true. A Runge-Kutta method should provide 
quicker access to the asymptotes than the Euler method. 

It was suggested that further work is intended on our simultaneous replacements and accelerated 
successive replacements iterative methods for nonlinear systems (ref. 3) and on another method which 
should almost converge unconditionally and derives from nonlinear least-squares considerations. 
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SOLUTION OF LINEAR SYSTEMS OF VERY LARGE SIZE 


S. W. Henriksen 

Raytheon Co. /Autometric Operation 

The subject of how to solve linear systems of very large size (high rank) is too large for a short 
discussion. A particular method of solution having outstanding peculiarities that are theoretically 
important will therefore be discussed. Idle method, invented by Michael Creusen (formerly of Auto- 
metric/Raytheon) is characterized by the following unusual features. 

(1) The customary system to be solved has the form 

Y _ A X m < n 

n x 1 ” nxm m x 1 


Creusen considers instead the system 



B Z 

nxm m x 1 


m > n 


which is made equivalent to the former. 


(2) If we write for the variances: 






where 


W = BZ 
-WC = AZ 
Z =Z + AZ 


and 


then 


[i] ^[iW 



AZ =)B ' 1 (I-B 1 C 1 )f(B 0 B 1 )Z 
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Now we find a unique value for C by minimizing, not as usual the scatter in AY but instead the scat- 
ter in &. This gives 


C = [ 0 I] 



so that 


AZ = 


- 7 -- 7' 1 
zw ww 


w 


(3) It can be shown that the Creu sen procedure always converges to a solution, even when A is 
singular or near singular. 

( 4 ) Writing the problem so that m >n allows great flexibility in programing the equations by 
minimizing storage space and by allowing sequential sdlution in steps of any size. 

The procedure has, on the other hand, certain disadvantages. These have not been fully worked 
out yet. 
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HEURISTIC INTEGRATION PROGRAM 


I. E. Perlin and L. J. Gallaher 
Georgialnstitute of Technology 

This discussion reports on a heuristic program for the numerical integration of systems of ordinary 
differential equations. 

The computer program was designed to incorporate the following features: ( 1 ) to meet a user 
specified accuracy; (2) to be problem independent (i.e., any number or degree of coupled differential 
equations); (3) to be self-optimizingwith respect to step size, order, and integration methods used; 
and (4) to exhibit learning so as to use a past-performance history to determine methods and orders to 
be used. 

The integration methods used are: (1) the method of Adams, Bashforth, and Moulton; (2) Butch- 
er’s formulas for the Stetter-Gragg-Butchermethod; (3) Cowell’s method of nth-order differences; and 
(4) Shank’s formulas for the Runge-Kutta method. 

The programs described here are written in both single ( 1 1 -decimal place) and double (22-decimal 
place) precision Algol for the B 5500. An executive procedure acts in an administrative and bookkeep- 
ing capacity for the various integration methods. This executive procedure keeps the performance 
histories according to problem types, determines performance effectiveness of the methods and orders, 
and chooses those to be used. 

Several kinds of problems and a large range of accuracies were used to exercise the program with 
the following results: 

(1) All methods perform well with no method or order always being exceptionally superior to 
The others 

(2) For different kinds of problems and accuracies different methods and orders prove more 
effective 

T (3) Learning takes place in a satisfactory manner 

(4) The program makes an effective general library procedure for integrating systems of ordinary 
differential equations. 
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NUMERICAL TECHNIQUES IN ORBITAL LIFETIME PREDICTIONS 


Edwin F. Fleischman, Jr. 

Marshall Space Flight Center, NASA 

Since the necessity of integrating over long periods of time, when one is calculating orbital life- 
time, the classical methods of trajectory calculations cannot be used. The rate of change of apogee 
and perigee is a slowly varying quantity, when compared to the change in some sort of position vec- 
tor. Hence, the rate of change of apogee and perigee can be integrated for long periods of time 
w hi le keeping error down to a minimum. 

A Simpson’s rule of integration is used to calculate values of A and P, which are orbit and 
vehicle dependent. These values are then integrated in time by a fourth-order Runga-Kutta scheme. 

Accuracy is within the state of art of the environmental modeling, and orbital ephemerides as 
long as 10 000 days can be calculated in a matter of minutes. 
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